Vamos a analizar de forma descriptiva algunas serie de tiempo. Empezaremos por la serie de desempleo de los Estados Unidos que viene en el paquete TSstudio.
library(TSstudio)
data(USUnRate)
ts_info(USUnRate)
## The USUnRate series is a ts object with 1 variable and 864 observations
## Frequency: 12
## Start time: 1948 1
## End time: 2019 12
class(USUnRate)
## [1] "ts"
plot(USUnRate,main = "US Monthly Unemployment Rate",ylab="Unemployment Rate (%)")
La serie es medida mensual, es decir, presenta una frecuencia de 12. Qué características podemos observar? * Tendencia? * Estacionalidad? * Cíclos? * Varianza marginal no constante? Vamos a seleccionar un periodo de tiempo mas corto
unemployment <- window(USUnRate, start = c(1990,1))
ts_plot(unemployment,
title = "US Monthly Unemployment Rate",
Ytitle = "Unemployment Rate (%)",
Xtitle = "Year",
Xgrid = TRUE,
Ygrid = TRUE)
Note que aquí podemos ver varias características: Estacionalidad(no es tan evidente), tres periodos de ciclos, el primero de 1990 a 2000,el segundo de 2000 a 2007, y el tercero de 2007 a 2019. No parece tener una heterocedasticidad marginal.
Veamos ahora la tasa de desempleo de Colombia. Hay que hacerle un ajuste a la base de datos porque está en orden descendente en el tiempo.
library(readxl)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
DesempleoyEmpleo <- read_excel("DesempleoyEmpleo.xlsx", range="A6:C269")
str(DesempleoyEmpleo)
## tibble [263 × 3] (S3: tbl_df/tbl/data.frame)
## $ Fecha : chr [1:263] "2022-11" "2022-10" "2022-09" "2022-08" ...
## $ Tasa_de_ocupacion: num [1:263] 57.4 57.7 57.2 56.7 56.5 ...
## $ Tasa_de_desempleo: num [1:263] 9.5 9.72 10.75 10.63 10.99 ...
DesempleoyEmpleo_1=DesempleoyEmpleo %>% map_df(rev)
tail(DesempleoyEmpleo)
## # A tibble: 6 × 3
## Fecha Tasa_de_ocupacion Tasa_de_desempleo
## <chr> <dbl> <dbl>
## 1 2001-06 56.2 15.3
## 2 2001-05 56.2 14.0
## 3 2001-04 55.8 14.5
## 4 2001-03 57.6 15.8
## 5 2001-02 56.9 17.4
## 6 2001-01 57.6 16.6
head(DesempleoyEmpleo_1)
## # A tibble: 6 × 3
## Fecha Tasa_de_ocupacion Tasa_de_desempleo
## <chr> <dbl> <dbl>
## 1 2001-01 57.6 16.6
## 2 2001-02 56.9 17.4
## 3 2001-03 57.6 15.8
## 4 2001-04 55.8 14.5
## 5 2001-05 56.2 14.0
## 6 2001-06 56.2 15.3
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(xts)
##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## ###############################################################################
##
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
##
## first, last
Fechas=as.yearmon(DesempleoyEmpleo_1$Fecha)
Desempleo_Col_xts=xts(x = DesempleoyEmpleo_1$Tasa_de_desempleo,frequency = 12,order.by = Fechas)
ts_info(Desempleo_Col_xts)
## The Desempleo_Col_xts series is a xts object with 1 variable and 263 observations
## Frequency: monthly
## Start time: Jan 2001
## End time: Nov 2022
plot(Desempleo_Col_xts)
ts_plot(Desempleo_Col_xts,
title = "Tasa de Desemplo Mensual Colombia",
Ytitle = "Tasa de Desempleo(%)",
Xtitle = "Año",
Xgrid = TRUE,
Ygrid = TRUE)
Qué características presenta esta serie?
Note que también se puede crear un objeto tsibble y usar el paquete feast, fable y fabletools(tsibble) y timetk(tibble), lo cuales funcionan con el paquete tidyverse.
require(feasts)
## Loading required package: feasts
## Loading required package: fabletools
require(fable)
## Loading required package: fable
require(timetk)
## Loading required package: timetk
require(tsibble)
## Loading required package: tsibble
##
## Attaching package: 'tsibble'
## The following object is masked from 'package:zoo':
##
## index
## The following object is masked from 'package:lubridate':
##
## interval
## The following objects are masked from 'package:base':
##
## intersect, setdiff, union
require(lubridate)
###Creación objeto tssible a partir de un objeto tibble
df_desempleo=data.frame(Desempleo=DesempleoyEmpleo_1$Tasa_de_desempleo,Fecha=DesempleoyEmpleo_1$Fecha)
tbl_desempleo=tibble(df_desempleo)
tbl_desempleo_format_fecha=tbl_desempleo
tbl_desempleo_format_fecha$Fecha=yearmonth(tbl_desempleo_format_fecha$Fecha)
###El tipo de fechas debe ser alguno que reconozca tsibble
tsbl_desempleo=as_tsibble(tbl_desempleo_format_fecha,index=Fecha) ####La fecha en tsibble es importante
##Gráfica de tsibble
autoplot(tsbl_desempleo,Desempleo)+labs(tittle="Serie de Desempleo Colombia Mensual",y="Tasa de Desempleo")
tbl_desempleo$Fecha<-as.Date(zoo::as.yearmon(tbl_desempleo$Fecha))
tbl_desempleo
## # A tibble: 263 × 2
## Desempleo Fecha
## <dbl> <date>
## 1 16.6 2001-01-01
## 2 17.4 2001-02-01
## 3 15.8 2001-03-01
## 4 14.5 2001-04-01
## 5 14.0 2001-05-01
## 6 15.3 2001-06-01
## 7 15.2 2001-07-01
## 8 14.9 2001-08-01
## 9 14.1 2001-09-01
## 10 14.5 2001-10-01
## # ℹ 253 more rows
###Gráfica timetk
tbl_desempleo%>%plot_time_series(.value=Desempleo,.date_var=Fecha)
Vamos a ver la forma de estimar la tendencia y/o eliminarla.
library(astsa)
library(TSstudio)
data(chicken)
ts_info(chicken)
## The chicken series is a ts object with 1 variable and 180 observations
## Frequency: 12
## Start time: 2001 8
## End time: 2016 7
plot(chicken,main="Precio Mensual de la Libra de Pollo en Estados Unidos", ylab="Precio en Centavos de Dólar")
#ts_plot(chicken)
Al parecer la serie de precios mensuales del pollo presenta una tendencia creciente al parecer lineal, es decir \[y_t=\mu_t+a_t\] o mas específicamente
\[y_t=\beta_0+\beta_1 t +a_t\]
Un primer intento de ver como sería la tendencia sería a través del uso del filtro de promedios móviles.
Vamos a usar la descomposición por medio de filtros de promedios móviles
chicken_decompo=decompose(chicken)
plot(chicken_decompo)
chicken_decompo$trend
## Jan Feb Mar Apr May Jun Jul
## 2001
## 2002 NA 63.91958 63.75667 63.54250 63.32875 63.15500 63.05458
## 2003 63.64792 63.96792 64.36708 64.81542 65.31667 65.89792 66.51458
## 2004 71.92417 72.96750 73.82292 74.50458 75.07875 75.53583 75.88917
## 2005 75.35458 74.87458 74.52875 74.33750 74.18583 74.00208 73.75333
## 2006 71.02208 70.63875 70.27000 69.88542 69.53458 69.30333 69.28708
## 2007 73.70708 74.62875 75.53333 76.40667 77.19292 77.87083 78.43000
## 2008 80.89500 81.48792 82.07125 82.68125 83.38750 84.19292 85.03333
## 2009 87.24208 87.18625 86.97083 86.62875 86.23333 85.83042 85.45208
## 2010 84.68500 84.69750 84.85958 85.14083 85.44125 85.71333 85.92833
## 2011 86.39167 86.38500 86.45042 86.59625 86.84792 87.19125 87.60042
## 2012 91.02625 91.62125 92.18625 92.74917 93.34292 93.97792 94.66958
## 2013 99.52667 100.49167 101.40917 102.23292 102.95292 103.56333 104.05833
## 2014 106.44125 106.96125 107.52875 108.20167 108.95417 109.73583 110.53667
## 2015 114.26083 114.50958 114.68250 114.76000 114.75208 114.69875 114.60333
## 2016 113.03167 NA NA NA NA NA NA
## Aug Sep Oct Nov Dec
## 2001 NA NA NA NA NA
## 2002 63.03542 63.09125 63.18125 63.27625 63.41958
## 2003 67.17167 67.90875 68.76083 69.72792 70.79583
## 2004 76.14000 76.26292 76.26458 76.13750 75.82708
## 2005 73.41375 72.99042 72.48750 71.95000 71.45333
## 2006 69.53958 70.06750 70.84500 71.77125 72.74708
## 2007 78.90083 79.32750 79.69417 80.02250 80.39333
## 2008 85.76458 86.26667 86.59333 86.87833 87.12667
## 2009 85.13500 84.92125 84.84500 84.81958 84.75667
## 2010 86.08375 86.24417 86.37750 86.42792 86.42208
## 2011 88.07750 88.61125 89.17625 89.77500 90.40333
## 2012 95.41000 96.14625 96.89542 97.70167 98.58000
## 2013 104.45875 104.79708 105.15125 105.54292 105.96083
## 2014 111.32708 112.08917 112.78208 113.39792 113.91000
## 2015 114.46792 114.28542 114.03375 113.72917 113.39000
## 2016
Simule unos datos III, crear un objeto con periodicidad s=12, y extraer la descomposición, qué puede observar?
y=ts(rnorm(1000,0,1),start=c(2000,01),frequency = 12)
plot(decompose(y))
Se puede calcular manualmente los filtros
chicken_decompo$trend
## Jan Feb Mar Apr May Jun Jul
## 2001
## 2002 NA 63.91958 63.75667 63.54250 63.32875 63.15500 63.05458
## 2003 63.64792 63.96792 64.36708 64.81542 65.31667 65.89792 66.51458
## 2004 71.92417 72.96750 73.82292 74.50458 75.07875 75.53583 75.88917
## 2005 75.35458 74.87458 74.52875 74.33750 74.18583 74.00208 73.75333
## 2006 71.02208 70.63875 70.27000 69.88542 69.53458 69.30333 69.28708
## 2007 73.70708 74.62875 75.53333 76.40667 77.19292 77.87083 78.43000
## 2008 80.89500 81.48792 82.07125 82.68125 83.38750 84.19292 85.03333
## 2009 87.24208 87.18625 86.97083 86.62875 86.23333 85.83042 85.45208
## 2010 84.68500 84.69750 84.85958 85.14083 85.44125 85.71333 85.92833
## 2011 86.39167 86.38500 86.45042 86.59625 86.84792 87.19125 87.60042
## 2012 91.02625 91.62125 92.18625 92.74917 93.34292 93.97792 94.66958
## 2013 99.52667 100.49167 101.40917 102.23292 102.95292 103.56333 104.05833
## 2014 106.44125 106.96125 107.52875 108.20167 108.95417 109.73583 110.53667
## 2015 114.26083 114.50958 114.68250 114.76000 114.75208 114.69875 114.60333
## 2016 113.03167 NA NA NA NA NA NA
## Aug Sep Oct Nov Dec
## 2001 NA NA NA NA NA
## 2002 63.03542 63.09125 63.18125 63.27625 63.41958
## 2003 67.17167 67.90875 68.76083 69.72792 70.79583
## 2004 76.14000 76.26292 76.26458 76.13750 75.82708
## 2005 73.41375 72.99042 72.48750 71.95000 71.45333
## 2006 69.53958 70.06750 70.84500 71.77125 72.74708
## 2007 78.90083 79.32750 79.69417 80.02250 80.39333
## 2008 85.76458 86.26667 86.59333 86.87833 87.12667
## 2009 85.13500 84.92125 84.84500 84.81958 84.75667
## 2010 86.08375 86.24417 86.37750 86.42792 86.42208
## 2011 88.07750 88.61125 89.17625 89.77500 90.40333
## 2012 95.41000 96.14625 96.89542 97.70167 98.58000
## 2013 104.45875 104.79708 105.15125 105.54292 105.96083
## 2014 111.32708 112.08917 112.78208 113.39792 113.91000
## 2015 114.46792 114.28542 114.03375 113.72917 113.39000
## 2016
as.numeric(chicken_decompo$trend)
## [1] NA NA NA NA NA NA 63.91958
## [8] 63.75667 63.54250 63.32875 63.15500 63.05458 63.03542 63.09125
## [15] 63.18125 63.27625 63.41958 63.64792 63.96792 64.36708 64.81542
## [22] 65.31667 65.89792 66.51458 67.17167 67.90875 68.76083 69.72792
## [29] 70.79583 71.92417 72.96750 73.82292 74.50458 75.07875 75.53583
## [36] 75.88917 76.14000 76.26292 76.26458 76.13750 75.82708 75.35458
## [43] 74.87458 74.52875 74.33750 74.18583 74.00208 73.75333 73.41375
## [50] 72.99042 72.48750 71.95000 71.45333 71.02208 70.63875 70.27000
## [57] 69.88542 69.53458 69.30333 69.28708 69.53958 70.06750 70.84500
## [64] 71.77125 72.74708 73.70708 74.62875 75.53333 76.40667 77.19292
## [71] 77.87083 78.43000 78.90083 79.32750 79.69417 80.02250 80.39333
## [78] 80.89500 81.48792 82.07125 82.68125 83.38750 84.19292 85.03333
## [85] 85.76458 86.26667 86.59333 86.87833 87.12667 87.24208 87.18625
## [92] 86.97083 86.62875 86.23333 85.83042 85.45208 85.13500 84.92125
## [99] 84.84500 84.81958 84.75667 84.68500 84.69750 84.85958 85.14083
## [106] 85.44125 85.71333 85.92833 86.08375 86.24417 86.37750 86.42792
## [113] 86.42208 86.39167 86.38500 86.45042 86.59625 86.84792 87.19125
## [120] 87.60042 88.07750 88.61125 89.17625 89.77500 90.40333 91.02625
## [127] 91.62125 92.18625 92.74917 93.34292 93.97792 94.66958 95.41000
## [134] 96.14625 96.89542 97.70167 98.58000 99.52667 100.49167 101.40917
## [141] 102.23292 102.95292 103.56333 104.05833 104.45875 104.79708 105.15125
## [148] 105.54292 105.96083 106.44125 106.96125 107.52875 108.20167 108.95417
## [155] 109.73583 110.53667 111.32708 112.08917 112.78208 113.39792 113.91000
## [162] 114.26083 114.50958 114.68250 114.76000 114.75208 114.69875 114.60333
## [169] 114.46792 114.28542 114.03375 113.72917 113.39000 113.03167 NA
## [176] NA NA NA NA NA
m_filter=12
filter_1=stats::filter(chicken, filter = rep(1/m_filter, m_filter), sides = 2)
filter_2=forecast::ma(chicken,m_filter)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
filter_3=forecast::ma(chicken,m_filter,centre = T)
Tarea: Estimar yEliminar la tendencia por un filtro de promedio móviles, y a la serie resultante extraerle la estacionalidad a partir de otro filtro de promedios móviles.(note que para este caso lo que hay que hacer es agrupar por meses a través de todos los años y calcular el promedio.)
Para mas detalles https://otexts.com/fpp3/moving-averages.html.
summary(fit <- lm(chicken~time(chicken), na.action=NULL))
##
## Call:
## lm(formula = chicken ~ time(chicken), na.action = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7411 -3.4730 0.8251 2.7738 11.5804
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.131e+03 1.624e+02 -43.91 <2e-16 ***
## time(chicken) 3.592e+00 8.084e-02 44.43 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.696 on 178 degrees of freedom
## Multiple R-squared: 0.9173, Adjusted R-squared: 0.9168
## F-statistic: 1974 on 1 and 178 DF, p-value: < 2.2e-16
plot(chicken, ylab="centavos por libra")
abline(fit,col = "red") # Se añade la recta ajusta
###Eliminamos la tendencia con la predicción la recta
ElimiTendchick=chicken-predict(fit)
plot(ElimiTendchick,main="Serie Chicken Sin tendencia")
acf(ElimiTendchick,lag.max =179 )
library(tidyverse)
library(lubridate)
library(timetk)
library(tsibble)
# Se configura para gráficas plotly(# FALSE retorna ggplots y no plotly)
interactive <- FALSE
indice_chicken=as.Date(as.yearmon(tk_index(chicken)))
## Otra forma de extraer el indice estimetk::tk_index(chicken)
df_chicken=data.frame(Fecha=indice_chicken,Pollo=as.matrix(chicken))
str(df_chicken)
## 'data.frame': 180 obs. of 2 variables:
## $ Fecha: Date, format: "2001-08-01" "2001-09-01" ...
## $ Pollo: num 65.6 66.5 65.7 64.3 63.2 ...
tibble_chicken=tibble(df_chicken)
duplicates(tibble_chicken, key = NULL, index=Fecha) ##Mirar si hay registros duplicados
## # A tibble: 0 × 2
## # ℹ 2 variables: Fecha <date>, Pollo <dbl>
tibble_chicken_fechas_correct=tibble_chicken
tibble_chicken_fechas_correct$Fecha=yearmonth(tibble_chicken_fechas_correct$Fecha)
print(duplicates(tibble_chicken_fechas_correct, key = NULL, index=Fecha))
## # A tibble: 0 × 2
## # ℹ 2 variables: Fecha <mth>, Pollo <dbl>
tsibble_chicken=tsibble(tibble_chicken_fechas_correct,index=Fecha)
tsibble_chicken
## # A tsibble: 180 x 2 [1M]
## Fecha Pollo
## <mth> <dbl>
## 1 2001 Aug 65.6
## 2 2001 Sep 66.5
## 3 2001 Oct 65.7
## 4 2001 Nov 64.3
## 5 2001 Dec 63.2
## 6 2002 Jan 62.9
## 7 2002 Feb 62.9
## 8 2002 Mar 62.7
## 9 2002 Apr 62.5
## 10 2002 May 63.4
## # ℹ 170 more rows
Note que se ajusta una tendencia que se ve que no es lineal.
tibble_chicken%>%timetk::plot_time_series(Fecha, Pollo,
.interactive = TRUE,
.plotly_slider = TRUE)
###Usa Loess para hacer el ajuste de la tendencia, es decir usar smooth_vec() como versión simplificada de stats::loess()
tibble_chicken%>%mutate(Pollo_ajus=smooth_vec(Pollo,span = 0.35, degree = 2))
## # A tibble: 180 × 3
## Fecha Pollo Pollo_ajus
## <date> <dbl> <dbl>
## 1 2001-08-01 65.6 62.8
## 2 2001-09-01 66.5 62.8
## 3 2001-10-01 65.7 62.9
## 4 2001-11-01 64.3 63.0
## 5 2001-12-01 63.2 63.2
## 6 2002-01-01 62.9 63.3
## 7 2002-02-01 62.9 63.5
## 8 2002-03-01 62.7 63.6
## 9 2002-04-01 62.5 63.8
## 10 2002-05-01 63.4 64.0
## # ℹ 170 more rows
tibble_chicken%>%mutate(Pollo_ajus=smooth_vec(Pollo,span = 0.75, degree = 2))%>%
ggplot(aes(Fecha, Pollo)) +
geom_line() +
geom_line(aes(y = Pollo_ajus), color = "red")
STL son las iniciales de “Seasonal and Trend decomposition using Loess”,el cual fue desarrollado por R. B. Cleveland et al. (1990).
library(feasts)
library(fable)
tsibble_chicken<-as_tsibble(chicken)
str(tsibble_chicken)
## tbl_ts [180 × 2] (S3: tbl_ts/tbl_df/tbl/data.frame)
## $ index: mth [1:180] 2001 Aug, 2001 Sep, 2001 Oct, 2001 Nov, 2001 Dec, 2002 Jan...
## $ value: num [1:180] 65.6 66.5 65.7 64.3 63.2 ...
## - attr(*, "key")= tibble [1 × 1] (S3: tbl_df/tbl/data.frame)
## ..$ .rows: list<int> [1:1]
## .. ..$ : int [1:180] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..@ ptype: int(0)
## - attr(*, "index")= chr "index"
## ..- attr(*, "ordered")= logi TRUE
## - attr(*, "index2")= chr "index"
## - attr(*, "interval")= interval [1:1] 1M
## ..@ .regular: logi TRUE
tsibble_chicken %>%
model(
STL(value ~ trend() +
season(window = "periodic"),
robust = TRUE)) %>%
components() %>%
autoplot()
Note que en ambos casos se obliga a extraer un componente estacional, sin embargo puede que está componente en verdad no exista, por eso se debe verificar que en efecto hay.
set.seed(154)
w = rnorm(200); x = cumsum(w)
wd = w +.2; xd = cumsum(wd)
plot.ts(xd, ylim=c(-5,55), main="Caminata Aletoria", ylab='')
lines(x, col=4); abline(h=0, col=4, lty=2); abline(a=0, b=.2, lty=2)
La diferencia ordinaria de orden 1 es , \[\nabla^1 Y_t=(1-B)^1 Y_t=Y_t-Y_{t-1}\]
dx=diff(x)
plot.ts(dx, ylim=c(-5,5), main="Serie Diferenciada", ylab='')
####Con drift
dxd=diff(xd)
plot.ts(dxd, ylim=c(-5,5), main="Serie con drift Diferenciada", ylab='')
par(mar = c(2,2,2,2))
fit = lm(chicken~time(chicken), na.action=NULL) # Regresión sobre el tiempo
par(mfrow=c(2,1))
plot(resid(fit), type="l", main="sin tendencia")
plot(diff(chicken), type="l", main="Primera Diferencia")
par(mar = c(3,2,3,2))
par(mfrow=c(3,1)) # plot ACFs
acf(chicken, 48, main="ACF Pollo")
acf(resid(fit), 48, main="ACF Sin tendencia")
acf(diff(chicken), 48, main="ACF Primera Diferencia")
En ocasiones la serie presenta varianza marginal no constante a lo largo del tiempo, lo cual hace necesario tener en cuenta tal característica. En este caso, se siguiere hacer una transformación de potencia para estabilizar la varianza. Esta familia de transformaciones se llaman transformaciones Box-Cox.
\[ f_{\lambda}(u_{t})= \begin{cases} \lambda^{-1}(u^{\lambda}_{t}-1), & \text{si $u_{t} \geq 0$, para $\lambda>0$,}\\ \ln(u_{t}), &\text{ si $u_{t}>0$, para $\lambda=0$}. \end{cases} \]
data("AirPassengers")
plot(AirPassengers)
#####Transformación Box-Cox
#library(FitAR)
library(forecast)
##
## Attaching package: 'forecast'
## The following object is masked from 'package:astsa':
##
## gas
## The following object is masked from 'package:fabletools':
##
## accuracy
forecast::BoxCox.lambda(AirPassengers, method ="loglik", lower = -1, upper = 3) ###Me entrega el valor de lambda
## [1] 0.2
##method="loglik"
#FitAR::BoxCox(AirPassengers)###Me entrega una gráfica
plot(forecast::BoxCox(AirPassengers,lambda=0.2))###
lAirPass=log(AirPassengers)
#x11()
par(mar = c(1,1,1,1))
par(mfrow=c(2,1))
plot(AirPassengers,main="Serie de Pasajeros sin Transformar")
plot(lAirPass,main="Series con Transformación BoxCox")
forecast::BoxCox.lambda(lAirPass, method ="guerrero", lower = -1, upper = 3)
## [1] -0.4190648
##Box-Cox con timetk
timetk::box_cox_vec(AirPassengers,lambda = 'auto',silent = F)
## box_cox_vec(): Using value for lambda: -0.294715585559316
## Jan Feb Mar Apr May Jun Jul Aug
## 1949 2.548484 2.561375 2.588408 2.582937 2.567506 2.593720 2.615089 2.615089
## 1950 2.555038 2.577299 2.603899 2.593720 2.575381 2.616631 2.646225 2.646225
## 1951 2.610379 2.618160 2.656279 2.636912 2.648795 2.656279 2.680103 2.680103
## 1952 2.647515 2.658701 2.673640 2.659900 2.662270 2.699009 2.709884 2.720049
## 1953 2.676904 2.676904 2.715050 2.714201 2.709007 2.720866 2.737089 2.742835
## 1954 2.685298 2.668053 2.714201 2.707236 2.713347 2.737089 2.762580 2.756933
## 1955 2.720049 2.712489 2.739270 2.740706 2.741419 2.770363 2.796341 2.787869
## 1956 2.751056 2.746317 2.771524 2.769193 2.772100 2.801088 2.818144 2.814820
## 1957 2.770363 2.761963 2.792420 2.788382 2.791921 2.821786 2.837892 2.838594
## 1958 2.784223 2.772100 2.795371 2.788382 2.795857 2.826872 2.846724 2.851232
## 1959 2.794394 2.785275 2.815241 2.810978 2.820985 2.840332 2.864126 2.867216
## 1960 2.819775 2.808794 2.820583 2.836477 2.840332 2.860370 2.883509 2.879580
## Sep Oct Nov Dec
## 1949 2.595457 2.563441 2.529834 2.561375
## 1950 2.629937 2.590196 2.552878 2.602242
## 1951 2.663443 2.635540 2.611963 2.640966
## 1952 2.690331 2.671428 2.648795 2.674735
## 1953 2.715895 2.692301 2.658701 2.682201
## 1954 2.733382 2.709007 2.684272 2.709007
## 1955 2.768604 2.744238 2.715895 2.747003
## 1956 2.791921 2.765020 2.742129 2.765020
## 1957 2.814399 2.787869 2.764414 2.782096
## 1958 2.814399 2.793903 2.767420 2.782631
## 1959 2.837187 2.815659 2.795371 2.814820
## 1960 2.852177 2.836477 2.808352 2.825716
Note que ahora usamos la misma función para verificar si en verdad la varianza fue estabilizada.
#FitAR::BoxCox(lAirPass)
forecast::BoxCox.lambda(lAirPass, method = "guerrero", lower = -1, upper = 2)
## [1] -0.419081
library(VGAM)
## Loading required package: stats4
## Loading required package: splines
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:VGAM':
##
## logit
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
VGAM::yeo.johnson(AirPassengers, lambda = 0)
## Jan Feb Mar Apr May Jun Jul Aug
## 1949 4.727388 4.779123 4.890349 4.867534 4.804021 4.912655 5.003946 5.003946
## 1950 4.753590 4.844187 4.955827 4.912655 4.836282 5.010635 5.141664 5.141664
## 1951 4.983607 5.017280 5.187386 5.099866 5.153292 5.187386 5.298317 5.298317
## 1952 5.147494 5.198497 5.267858 5.204007 5.214936 5.389072 5.442418 5.493061
## 1953 5.283204 5.283204 5.468060 5.463832 5.438079 5.497168 5.579730 5.609472
## 1954 5.323010 5.241747 5.463832 5.429346 5.459586 5.579730 5.713733 5.683580
## 1955 5.493061 5.455321 5.590987 5.598422 5.602119 5.755742 5.899897 5.852202
## 1956 5.652489 5.627621 5.762051 5.749393 5.765191 5.926926 6.025866 6.006353
## 1957 5.755742 5.710427 5.877736 5.855072 5.874931 6.047372 6.144186 6.148468
## 1958 5.831882 5.765191 5.894403 5.855072 5.897154 6.077642 6.198479 6.226537
## 1959 5.888878 5.837730 6.008813 5.983936 6.042633 6.159095 6.308098 6.327937
## 1960 6.035481 5.971262 6.040255 6.135565 6.159095 6.284134 6.434547 6.408529
## Sep Oct Nov Dec
## 1949 4.919981 4.787492 4.653960 4.779123
## 1950 5.068904 4.897840 4.744932 4.948760
## 1951 5.220356 5.093750 4.990433 5.117994
## 1952 5.347108 5.257495 5.153292 5.273000
## 1953 5.472271 5.356586 5.198497 5.308268
## 1954 5.560682 5.438079 5.318120 5.438079
## 1955 5.746203 5.616771 5.472271 5.631212
## 1956 5.874931 5.726848 5.605802 5.726848
## 1957 6.003887 5.852202 5.723585 5.820083
## 1958 6.003887 5.886104 5.739793 5.823046
## 1959 6.139885 6.011267 5.894403 6.006353
## 1960 6.232448 6.135565 5.968708 6.070738
car::yjPower(AirPassengers,lambda=0)
## [1] 4.727388 4.779123 4.890349 4.867534 4.804021 4.912655 5.003946 5.003946
## [9] 4.919981 4.787492 4.653960 4.779123 4.753590 4.844187 4.955827 4.912655
## [17] 4.836282 5.010635 5.141664 5.141664 5.068904 4.897840 4.744932 4.948760
## [25] 4.983607 5.017280 5.187386 5.099866 5.153292 5.187386 5.298317 5.298317
## [33] 5.220356 5.093750 4.990433 5.117994 5.147494 5.198497 5.267858 5.204007
## [41] 5.214936 5.389072 5.442418 5.493061 5.347108 5.257495 5.153292 5.273000
## [49] 5.283204 5.283204 5.468060 5.463832 5.438079 5.497168 5.579730 5.609472
## [57] 5.472271 5.356586 5.198497 5.308268 5.323010 5.241747 5.463832 5.429346
## [65] 5.459586 5.579730 5.713733 5.683580 5.560682 5.438079 5.318120 5.438079
## [73] 5.493061 5.455321 5.590987 5.598422 5.602119 5.755742 5.899897 5.852202
## [81] 5.746203 5.616771 5.472271 5.631212 5.652489 5.627621 5.762051 5.749393
## [89] 5.765191 5.926926 6.025866 6.006353 5.874931 5.726848 5.605802 5.726848
## [97] 5.755742 5.710427 5.877736 5.855072 5.874931 6.047372 6.144186 6.148468
## [105] 6.003887 5.852202 5.723585 5.820083 5.831882 5.765191 5.894403 5.855072
## [113] 5.897154 6.077642 6.198479 6.226537 6.003887 5.886104 5.739793 5.823046
## [121] 5.888878 5.837730 6.008813 5.983936 6.042633 6.159095 6.308098 6.327937
## [129] 6.139885 6.011267 5.894403 6.006353 6.035481 5.971262 6.040255 6.135565
## [137] 6.159095 6.284134 6.434547 6.408529 6.232448 6.135565 5.968708 6.070738
Ejercicio: Use la serie varve del paquete astsa para chequear si es necesario hacer transformación Box-Cox.
Vamos a correr lo mismo pero en Python para la transformación de BoxCox
#```{r configuración, include=FALSE} #library(knitr) #library(reticulate) ###Nos ubicamos en el ambiente dentro de la terminal y damos which python y copiamos la ruta “/opt/anaconda3/envs/Python38andR/bin/python” ##“/Users/macbook/opt/anaconda3/envs/Ambiente38/bin/python” #/Users/sergiocalderonunal/opt/anaconda3/envs/Ambiente3104/bin/python #use_python(“/Users/macbook/opt/anaconda3/envs/Ambiente38/bin/python”) #use_python(“/Users/sergiocalderonunal/opt/anaconda3/envs/Ambiente3104/bin/python”) #use_virtualenv(“~/Python38andR”) #py_config() #pd <- import(“pandas”) #arreglo<-pd\(array(c(1, 2, 3)) #print(arreglo) #arreglo\)shape
#```
#```{python inicio} #import sys #print(sys.path) #import pandas as pd #import matplotlib.pylab as plt #data = pd.read_csv(‘AirPassengers.csv’) #print(data) #print(‘Data Types:’) #print(data.dtypes) #con=data[‘Month’] #data[‘Month’]=pd.to_datetime(data[‘Month’]) #data #pasajeros=data.set_index(‘Month’) #check datatype of index
#convert to time series: #ts = pasajeros[‘NPassengers’] #ts.head(10)
####Graficar la Serie##### #plt.plot(ts) #plt.title(‘AirPassengers’) #plt.show()
#```
#{python inicio_1} #print(r.AirPassengers) #print(r.lAirPass) #
#```{r uso de objeto creado en python}
#py\(data #library(ggplot2) #ggplot2::ggplot(data = py\)data,aes(x=Month,y=NPassengers) #)+geom_line()
#```
Vamos a hacer gráficos de dispersión para chequear que tipos de relaciones hay entre los retardos de la variable interés. Vamos a trabajar con algunas series, por ejemplo: * Indice ambiental mensual (soi)(Southern Oscillation Index), el cual mide los cambios en la presión del aire, relacionados con las temperaturas de la superficie del mar en el Océano Pacífico central. * la serie rec (reclutamiento asociada al soi), número de nuevos peces. * Consumo mensual de gas natural en EE. UU.(USgas) medido en Billones de pies cúbicos. Esto permite chequear si hay posibles relaciones no-lineales.
library(astsa)
data("soi")
ts_info(soi)
## The soi series is a ts object with 1 variable and 453 observations
## Frequency: 12
## Start time: 1950 1
## End time: 1987 9
?soi
par(mar = c(2,2,2,2))
plot(soi, main="Indice soi")
par(mar = c(3,2,3,2))
astsa::lag1.plot(soi, 12,corr=F) ###El 12 indica cuantos retardos y_t-k contra y_t
###Hacer la gráfica con x11()
En el gráfico de dispersión podemos ver que se muestra un ajuste no paramétrico, de la posible relación entre las variables al igual que una estimación de la autocorrelación entre \(s_t\) y \(s_{t-h}\). Vemos que varias de las relaciones exploradas parecen ser lineales. Uno puede observar que existen relaciones lineales positivas en los rezagos \(h = 1, 2, 11, 12\), mientras que negativas en los rezagos \(h=6,7\), las demás parecen ser no significativas o no lineales.
#pdf('/Users/sergiocalderon/Documents/Documentos - iMac de Sergio/Documentos iMac Sergio/Notas de Clase/Notas de clase/Notas de Clase Series de Tiempo Univariadas/Graficas/DispersionSoiRec.pdf',paper="USr")
par(mar = c(3,2,3,2))
lag2.plot(soi, rec, 8) #El 2 de lag2.plot es porque intervenienen dos serie de tiempo.
#dev.off()
lag2.plot(soi, rec, 8,corr=F)
Note también que el gráfico de dispersión del índice de nuevos peces con retardos del indice soi nos muestra también relaciones posiblemente lineales, al igual que no lineales. Ahora, podemos crear el gráfico de autocorrelación simple que nos permite estimar y graficar las autocorrelaciones para diferentes rezagos.
Note que con la función ts_lags de TSstudio podemos hacer una gráfica similar.
ts_lags(soi,lags=1:12)
Cuando el proceso es estacionario, o al menos no presenta tendencia, podemos usar el gráfico acf para explorar las posibles relaciones lineales a diferentes rezagos. En seguida mostramos la función de autocorrelación para el índice soi y la serie de nuevos peces.
par(mfrow=c(2,1))
par(mar = c(2.7,2,2.7,2))
acf(soi, 48, main="Southern Oscillation Index")
acf(rec, 48, main="Recruitment")
#dev.off()
ccf(rec, soi, 48, main="SOI vs Recruitment", ylab="CCF")
#Corr(rec_{t},soi_{t-h})
#Corr(soi_{t},rec_{t-h},)
index_soi_rec=as.Date(as.yearmon(tk_index(soi)))
df_soi_rec=data.frame(Fecha_soi_rec=index_soi_rec,soi=as.matrix(soi),rec=as.matrix(rec))
tibble_soi_rec=tibble(df_soi_rec)
tsibble::duplicates(tibble_soi_rec, key = NULL, index=Fecha_soi_rec) ##Mirar si hay registros duplicados
## # A tibble: 0 × 3
## # ℹ 3 variables: Fecha_soi_rec <date>, soi <dbl>, rec <dbl>
tibble_soi_rec%>%plot_acf_diagnostics(Fecha_soi_rec,soi,.ccf_vars = rec,.lags = 36)
En ambas gráficas podemos ver que se muestran periodicidades en las correlaciones que corresponden a valores separados por 12 unidades. También podemos ver que las observaciones con 12 meses o un año de diferencia están fuertemente correlacionadas positivamente, al igual que las observaciones en múltiplos como 24, 36, 48,\(\cdots\) Las observaciones separadas por seis meses están correlacionadas negativamente, lo que muestra que las excursiones positivas tienden a asociarse con las excursiones negativas a los seis meses ver Shumway2017 capítulo 1.
###AMI Del los libros H. Kantz and T. Schreiber: Nonlinear Time series Analysis (Cambridge university press) H. Abarbanel: Analysis of observed chaotic data (Springer, 1996) y NONLINEAR TIME SERIES ANALYSIS HOLGER KANTZ AND THOMAS SCHREIBER. Cambrige University Press 2003.
Ahora utilizaremos los paquetes nonlinearTseries y tseriesChaos para computar el average mutual information(AMI) o La información mutua promedio (AMI, la cual mide cuánto nos dice una variable aleatoria sobre otra, el cual se define como:
\[I(X;Y)=\sum_{i}\sum_{j}p(x_i,y_j)\log_2(\frac{p(x_i,y_j)}{p(x_i)p(y_j)}).\] En el contexto del análisis de series de tiempo, AMI ayuda a cuantificar la cantidad de conocimiento obtenido sobre el valor de \(X_{t+d}\) al observar \(X_t\). Equivalentemente, el AMI es una medida de qué tanto el conocimiento de \(X\) reduce la incertidumbre acerca de \(Y\). Esto implica que \(I(X,Y)=0\) si y sólo si \(X\) y \(Y\) son variables aletorias independientes. I(X; Y ) describe la información que la medición \(X_t\) en el tiempo \(t\) aporta a la medición \(X_{t+d}\) en el tiempo \(t + d\). Si se elige d como el valor alrededor del primer mínimo del AMI, entonces \(Y{t}\) e \(y_{t+d}\) son parcialmente pero no totalmente independientes.
Vamos a simular una serie de la forma \[x_t=\frac{x_{t-1}}{x_{t-12}^2+1}+\epsilon_t\] y a trabajar con la serie de linces Candienses del paquete astsa lynx.
library(nonlinearTseries)
##
## Attaching package: 'nonlinearTseries'
## The following object is masked from 'package:fabletools':
##
## estimate
## The following object is masked from 'package:grDevices':
##
## contourLines
library(tseriesChaos)
et=rnorm(1100,0,1)
xt=rep(0,1100)
for(t in 13:1100)
{
xt[t]=(xt[t-1]-1)/(xt[t-12]^2+1)+et[t]
}
xtsimul=as.ts(xt[101:1100])
length(xtsimul)
## [1] 1000
plot(xtsimul)
acf(xtsimul)
par(mar = c(3,2,3,2))
astsa::lag1.plot(xtsimul, 12,corr = F)
nonlinearTseries::mutualInformation(xtsimul,lag.max = 100,n.partitions = 50,units = "Bits",do.plot = TRUE) #c("Nats", "Bits", "Bans")
## $time.lag
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## [19] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
## [37] 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [55] 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
## [73] 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
## [91] 90 91 92 93 94 95 96 97 98 99 100
##
## $mutual.information
## [1] 4.7211355 0.9086685 0.7868768 0.7146046 0.7234251 0.7393485 0.7257831
## [8] 0.7204692 0.7530212 0.7388565 0.6987539 0.7426234 0.8363724 0.7918314
## [15] 0.7495687 0.7384157 0.7227720 0.7306765 0.7185728 0.7490403 0.7888126
## [22] 0.7335948 0.7676440 0.7241823 0.7728012 0.7690717 0.7278799 0.7436484
## [29] 0.7654490 0.7561503 0.7324961 0.7820421 0.7413698 0.7871444 0.7796870
## [36] 0.7356812 0.7411850 0.7501907 0.7512091 0.7691925 0.7466313 0.7425748
## [43] 0.7152309 0.7567790 0.7658219 0.7516127 0.7371519 0.7288299 0.7622234
## [50] 0.7844343 0.7584129 0.7388515 0.7564102 0.7690676 0.7334833 0.7605943
## [57] 0.7808753 0.7568963 0.7571249 0.7135538 0.7925615 0.7584517 0.7768902
## [64] 0.7298138 0.7462323 0.7553250 0.7728154 0.7850680 0.7464423 0.7543733
## [71] 0.7465533 0.7813147 0.7620542 0.7528362 0.7663679 0.8111562 0.7237969
## [78] 0.7386509 0.7523882 0.7398758 0.7844187 0.7326143 0.7380267 0.7455960
## [85] 0.7691433 0.7309565 0.7479580 0.7314371 0.7583339 0.8090117 0.7704502
## [92] 0.7519965 0.7847387 0.7292170 0.7291928 0.7169975 0.7426216 0.7759557
## [99] 0.7898280 0.7688683 0.7700155
##
## $units
## [1] "Bits"
##
## $n.partitions
## [1] 50
##
## attr(,"class")
## [1] "mutualInf"
pacf(xtsimul)
tseriesChaos::mutual(xtsimul, partitions = 50, lag.max = 100, plot=TRUE)
plot(astsa::Lynx)
acf(astsa::Lynx)
par(mar = c(3,2,3,2))
astsa::lag1.plot(astsa::Lynx, 12)
tseriesChaos::mutual(astsa::Lynx, partitions = 50, lag.max = 10, plot=TRUE)
Podemos usar la función TSstudio::ts_heatmap para crear un mapa de calor. Este es un gráfico tridimensional, en donde en el eje y están los meses, y en el eje x están los años. Note en este caso que los meses de Diciembre,Enero, Febrero y Marzo presentan los valores mas oscuros, es decir los valores mas grandes a lo largo de los años, en contraste con los meses de Mayo a Septiembre que presentan colores mas claros. Este es un típico comportamiento de la presencia de un cíclo estacional en la serie. El flujo de color es horizontal.
TSstudio::USgas
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 2000 2510.5 2330.7 2050.6 1783.3 1632.9 1513.1 1525.6 1653.1 1475.0 1567.8
## 2001 2677.0 2309.5 2246.6 1807.2 1522.4 1444.4 1598.1 1669.2 1494.1 1649.1
## 2002 2487.6 2242.4 2258.4 1881.0 1611.5 1591.4 1748.4 1725.7 1542.2 1645.9
## 2003 2700.5 2500.3 2197.9 1743.5 1514.7 1368.4 1600.5 1651.6 1428.6 1553.2
## 2004 2675.8 2511.1 2100.9 1745.2 1573.0 1483.7 1584.9 1578.0 1482.2 1557.2
## 2005 2561.9 2243.0 2205.8 1724.9 1522.6 1534.1 1686.6 1695.1 1422.5 1428.2
## 2006 2165.3 2144.4 2126.4 1681.0 1526.3 1550.9 1758.7 1751.7 1462.1 1644.2
## 2007 2475.6 2567.0 2128.8 1810.1 1559.1 1555.2 1659.9 1896.1 1590.5 1627.8
## 2008 2734.0 2503.4 2278.2 1823.9 1576.4 1604.2 1708.6 1682.9 1460.9 1635.8
## 2009 2729.7 2332.5 2170.7 1741.3 1504.0 1527.8 1658.0 1736.5 1575.0 1666.5
## 2010 2809.8 2481.0 2142.9 1691.8 1617.3 1649.5 1825.8 1878.9 1637.5 1664.9
## 2011 2888.6 2452.4 2230.5 1825.0 1667.4 1657.3 1890.5 1891.8 1655.6 1744.5
## 2012 2756.2 2500.7 2127.8 1953.1 1873.8 1868.4 2069.8 2008.8 1807.2 1901.1
## 2013 2878.8 2567.2 2521.1 1967.5 1752.5 1742.9 1926.3 1927.4 1767.0 1866.8
## 2014 3204.1 2741.2 2557.9 1961.7 1810.2 1745.4 1881.0 1933.1 1809.3 1912.8
## 2015 3115.0 2925.2 2591.3 2007.9 1858.1 1899.9 2067.7 2052.7 1901.3 1987.3
## 2016 3091.7 2652.3 2356.3 2083.8 1965.8 2000.7 2186.6 2208.4 1947.8 1925.2
## 2017 2914.2 2340.6 2523.7 1932.5 1892.5 1910.9 2142.1 2094.3 1920.9 2032.0
## 2018 3335.0 2705.9 2792.6 2346.3 2050.9 2058.7 2344.6 2307.7 2151.5 2279.1
## 2019 3399.9 2999.2 2899.9 2201.1 2121.0 2115.2 2407.5 2437.2 2215.6 2472.3
## Nov Dec
## 2000 1908.5 2587.5
## 2001 1701.0 2120.2
## 2002 1913.6 2378.9
## 2003 1753.6 2263.7
## 2004 1782.8 2327.7
## 2005 1663.4 2326.4
## 2006 1765.4 2122.8
## 2007 1834.5 2399.2
## 2008 1868.9 2399.7
## 2009 1776.2 2491.9
## 2010 1973.3 2714.1
## 2011 2031.9 2541.9
## 2012 2167.8 2503.9
## 2013 2316.9 2920.8
## 2014 2357.5 2679.2
## 2015 2249.1 2588.2
## 2016 2159.4 2866.3
## 2017 2357.7 3084.5
## 2018 2709.9 2993.1
## 2019
ts_plot(USgas,
title = "US Monthly Natural Gas consumption",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year",
Xgrid = TRUE,
Ygrid = TRUE)
TSstudio::ts_heatmap(USgas,title = "Mapa de Calor - Consumo de Gas Natural en EEUU")
Voy enfocarme en un periodo de la serie de desempleo.
unemployment <- window(USUnRate, start = c(1990,1))
ts_plot(unemployment,
title = "US Monthly Unemployment Rate",
Ytitle = "Unemployment Rate (%)",
Xtitle = "Year",
Xgrid = TRUE,
Ygrid = TRUE)
ts_info(unemployment)
## The unemployment series is a ts object with 1 variable and 360 observations
## Frequency: 12
## Start time: 1990 1
## End time: 2019 12
ts_heatmap(unemployment,title = "Mapa de Calor - Tasa de Desempleo EEUU Subserie")
#ts_heatmap(USUnRate,title = "Mapa de Calor - Tasa de Desempleo EEUU")
En este ejemplo, el flujo de color de la Tasa de Desempleo es vertical, lo que indica el estado del ciclo. En este caso, las franjas verticales más claras representan el final de un ciclo y el comienzo del siguiente. Asimismo, las franjas verticales más oscuras representan los picos del ciclo.
Vale la pena decir que es necesario eliminar la tendencia de la serie para pasar a detectar la estacionalidad.
USgas_df <- data.frame(year = floor(time(USgas)), month = cycle(USgas),USgas = as.numeric(USgas))
USgas_df$month <- factor(month.abb[USgas_df$month], levels = month.abb)
library(dplyr)
USgas_summary <- USgas_df %>%group_by(month) %>%summarise(mean= mean(USgas),sd = sd(USgas))
USgas_summary
## # A tibble: 12 × 3
## month mean sd
## <fct> <dbl> <dbl>
## 1 Jan 2806. 309.
## 2 Feb 2502. 223.
## 3 Mar 2325. 241.
## 4 Apr 1886. 175.
## 5 May 1708. 194.
## 6 Jun 1691. 216.
## 7 Jul 1864. 261.
## 8 Aug 1889. 237.
## 9 Sep 1687. 242.
## 10 Oct 1788. 260.
## 11 Nov 2015. 284.
## 12 Dec 2543. 278.
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
plot_ly (data = USgas_summary, x = ~ month, y = ~ mean, type = "bar", name = "Mean") %>%
layout (title = "USgas - Monthly Average", yaxis =list(title = "Mean", range = c(1500, 2700)))
monthplot(USgas)
####También lo podemos hacer usando el paquete feasts y el objeto tsibble
USgas_tsbl=as_tsibble(USgas)
require(feasts)
USgas_tsbl=as_tsibble(USgas)
USgas_tsbl%>%gg_subseries(value)###Puede usar el argumento period=12 y da el mismo resultado, lo que significa es que se pueden agrupar las observaciones que están cada 12.
Note que basados en las estadísticas descriptivas podemos ver que las medias son distintas para algunos meses, incluso sin quitar la tendencia. A una misma conclusión llegamos basados en los gráficos. Estas son típicas características de que hay presente un ciclo estacional en la serie.
Pueden haber múltiples estacionalidades(por lo general ocurren en series de alta frecuencia: diaria, cada hora, cada media y asi sucesivamente.) o una única estacionalidad, o inclusive cíclos que no se ven con facilidad.
Vamos a trabajar la base datos relacionada con el sistema nacional de transmisión de electricidad del Reino Unido. Más específicamente con la demanda de energía de forma horaria.
library(UKgrid)
require(TSstudio)
require(timetk)
require(feasts)
require(tsibble)
require(plotly)
UKgrid_xts <- extract_grid(type = "xts",
columns = "ND",
aggregate = "hourly",
na.rm = TRUE)
#extract_grid solo funciona para el conjunto de datos UKgrid
ts_plot(UKgrid_xts,
title = "National Hourly Demand UK Grid",
Ytitle = "Megawatts",
Xtitle = "Year",
Xgrid = TRUE,
Ygrid = TRUE)
UKgrid_tstbl <- extract_grid(type = "tsibble",
columns = "ND",
aggregate = "hourly",
na.rm = TRUE)
UKgrid_tbl <-na.omit(as_tibble(UKgrid_tstbl))
Como indicamos anteriormente, la primera indicación de la posible existencia de múltiples patrones estacionales en la serie es una frecuencia alta, como por diaria, horaria y en minutos. En esos casos, hay más de una forma de establecer la frecuencia de la serie. Por ejemplo, si capturamos una serie de tiempo de frecuencia diaria, la frecuencia de la serie se puede configurar de la siguiente manera: * Diariamente (o 365), asumiendo que el ciclo más apropiado es un año completo. * Días de semana (o 7) siempre que la oscilación del día de la semana sea más dominante que la del ciclo de año completo.
Al utilizar estadísticas descriptivas con este tipo de series, tendrá sentido aplicar este método para cada frecuencia potencial de la serie (o al menos las principales) con el fin de examinar si existe una indicación del patrón estacional.
Por ejemplo, UKgrid es una serie de tiempo por horas, que la marca automáticamente como sospechosa de tener múltiples patrones estacionales. Potencialmente, como se mencionó anteriormente, la demanda horaria de electricidad podría tener tres patrones estacionales diferentes:
Usando el paquete lubridate crearemos las características.
library(xts)
UKgrid_df <- data.frame(time = zoo::index(UKgrid_xts), UKgrid=as.numeric(UKgrid_xts))
str(UKgrid_df)
## 'data.frame': 127296 obs. of 2 variables:
## $ time : POSIXct, format: "2005-04-01 00:00:00" "2005-04-01 01:00:00" ...
## $ UKgrid: num 65080 68207 69172 66769 64660 ...
Ahora crearemos características estacionales basados en los periodos que deseamos explorar, por ejemplo hora del día,o día de la semana, o mes del año.
library(lubridate)
UKgrid_df$hour <- hour(UKgrid_df$time)
UKgrid_df$weekday <- wday(UKgrid_df$time, label = TRUE, abbr = TRUE)
UKgrid_df$month <- factor(month.abb[month(UKgrid_df$time)], levels = month.abb)
head(UKgrid_df)
## time UKgrid hour weekday month
## 1 2005-04-01 00:00:00 65080 0 Fri Apr
## 2 2005-04-01 01:00:00 68207 1 Fri Apr
## 3 2005-04-01 02:00:00 69172 2 Fri Apr
## 4 2005-04-01 03:00:00 66769 3 Fri Apr
## 5 2005-04-01 04:00:00 64660 4 Fri Apr
## 6 2005-04-01 05:00:00 65209 5 Fri Apr
Vamos a empezar las exploraciones analizando el ciclo horario.
UKgrid_hourly <- UKgrid_df %>%
dplyr::group_by(hour) %>%
dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE), sd = sd(UKgrid, na.rm
= TRUE))
str(UKgrid_hourly)
## tibble [24 × 3] (S3: tbl_df/tbl/data.frame)
## $ hour: int [1:24] 0 1 2 3 4 5 6 7 8 9 ...
## $ mean: num [1:24] 55622 54343 52920 51546 50509 ...
## $ sd : num [1:24] 9259 9482 9485 9165 8557 ...
UKgrid_hourly
## # A tibble: 24 × 3
## hour mean sd
## <int> <dbl> <dbl>
## 1 0 55622. 9259.
## 2 1 54343. 9482.
## 3 2 52920. 9485.
## 4 3 51546. 9165.
## 5 4 50509. 8557.
## 6 5 51747. 8657.
## 7 6 59048. 10489.
## 8 7 68627. 13347.
## 9 8 73581. 13003.
## 10 9 76320. 12557.
## # ℹ 14 more rows
Vamos ahora a hacer la gráfica de la media y la desviación estándar con base en las horas. Note que esas gráficas en escalas diferentes, así que hay que usar un gráfico especial.
require(plotly)
plot_ly(UKgrid_hourly) %>%
add_lines(x = ~ hour, y = ~ mean, name = "Media") %>%
add_lines(x = ~ hour, y = ~ sd, name = "Desviación Estándar", yaxis =
"y2",
line = list(color = "red", dash = "dash", width = 3)) %>%
layout(
title = "La demanda nacional de electricidad - Promedio horario vs. Desviación Estándar",
yaxis = list(title = "Media"),
yaxis2 = list(overlaying = "y",
side = "right",
title = "Desviación Estándar"
),
xaxis = list(title="Hora del Día"),
legend = list(x = 0.05, y = 0.9),
margin = list(l = 50, r = 50)
)
Qué podemos destacar del gráfico y de las estadísticas descriptivas?
Hay poca demanda durante la noche (entre la medianoche y las 6 a.m.) y una alta demanda entre las horas de la mañana y la tarde.
Existe una fuerte correlación entre la demanda promedio y su desviación estándar.
La relativamente baja desviación estándar de la demanda promedio durante la noche podría indicar que existe un fuerte efecto subestacional durante esas horas además de la estacionalidad horaria. Esto debería tener sentido, ya que son horas de sueño normales y, por lo tanto, en promedio, la demanda es razonablemente la misma durante los días de semana.
Por otro lado, la alta desviación estándar a lo largo de las horas de alta demanda podría indicar que la demanda se distribuye de manera diferente en diferentes vistas de periodicidad (como día de la semana o mes del año).
Vamos a explorar el último punto, viendo la demanda a la madrugada(3 a.m) y empezando el día(9 a.m) con respecto al día de la semana.
UKgrid_weekday <- UKgrid_df %>%
dplyr::filter(hour == 3 | hour == 9) %>%
dplyr::group_by(hour, weekday) %>%
dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE),
sd = sd(UKgrid, na.rm = TRUE))
## `summarise()` has grouped output by 'hour'. You can override using the
## `.groups` argument.
UKgrid_weekday$hour <- factor(UKgrid_weekday$hour)
plot_ly(data = UKgrid_weekday, x = ~ weekday, y = ~ mean, type =
"bar",color = ~ hour) %>%
layout(title = "The Hourly Average Demand by Weekday",
yaxis = list(title = "Mean", range = c(30000, 75000)),
xaxis = list(title = "Weekday"))
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
En el gráfico de barras anterior podemos ver que la demanda de electricidad a las 3 a.m. es relativamente estable durante todos los días de la semana, con una ligera diferencia entre el promedio durante los días laborables y los días del fin de semana (alrededor de un 2% diferente). Por otro lado, existe una diferencia entre la demanda del día laborable y el fin de semana a las 9 a.m. (es decir, la demanda del lunes es en promedio un 28% superior a la del domingo). Como era de esperar, esos resultados se alinearon con nuestras expectativas anteriores.
Ahora podemos aprovechar esos conocimientos para examinar si existe un patrón estacional mensual en la serie. Ahora seleccionaremos las mismas horas (3 a.m. y 9 a.m.); sin embargo, esta vez agruparemos estos datos por mes (en lugar de días laborables):
UKgrid_month <- UKgrid_df %>%
dplyr::filter(hour == 3 | hour == 9) %>%
dplyr::group_by(hour, month) %>%
dplyr::summarise(mean = mean(UKgrid, na.rm = TRUE),
sd = sd(UKgrid, na.rm = TRUE))
## `summarise()` has grouped output by 'hour'. You can override using the
## `.groups` argument.
UKgrid_month$hour <- factor(UKgrid_month$hour)
plot_ly(data = UKgrid_month, x = ~ month, y = ~ mean, type = "bar",color =
~ hour) %>%
layout(title = "The Hourly Average Demand by Month",
yaxis = list(title = "Mean", range = c(30000, 75000)),
xaxis = list(title = "Month"))
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Podemos ver en el gráfico de barras del resumen de agregación mensual que, en promedio, la demanda durante la noche (3 a.m.) y la mañana (9 a.m.) varía a lo largo de los meses del año. Además, hay un cambio significativo en la demanda durante la noche en comparación con la agregación entre semana. La variación de la serie de mes a mes indica la existencia de estacionalidad mensual en la serie.
Otro enfoque para analizar patrones estacionales en datos de series de tiempo es trazar la distribución de las unidades de frecuencia mediante el uso de histogramas o gráficos de densidad. Esto nos permitirá examinar si cada unidad de frecuencia tiene una distribución única que puede distinguirla del resto de unidades. Para esto usaremos el paquete ggplot2. Vamos empezar con la serie USgas.
require(timetk)
UKgrid_tbl
## # A tibble: 127,282 × 2
## TIMESTAMP ND
## <dttm> <int>
## 1 2005-04-01 00:00:00 65080
## 2 2005-04-01 01:00:00 68207
## 3 2005-04-01 02:00:00 69172
## 4 2005-04-01 03:00:00 66769
## 5 2005-04-01 04:00:00 64660
## 6 2005-04-01 05:00:00 65209
## 7 2005-04-01 06:00:00 72376
## 8 2005-04-01 07:00:00 82679
## 9 2005-04-01 08:00:00 88706
## 10 2005-04-01 09:00:00 91108
## # ℹ 127,272 more rows
#UKgrid_tbl%>%plot_seasonal_diagnostics(.date_var = TIMESTAMP,.value = ND,.feature_set = c("hour","week","wday.lbl","month.lbl"),.geom="boxplot") ##Chequear si se desea un diagrama de caja o violín .geom
UKgrid_diff<-diff(UKgrid_tbl$ND)
which(is.na(UKgrid_diff))
## integer(0)
UKgrid_diff_sub=na.omit(UKgrid_diff)
str(UKgrid_diff_sub)
## int [1:127281] 3127 965 -2403 -2109 549 7167 10303 6027 2402 209 ...
UKgrid_tbl%>%mutate(diff_ND=ND-dplyr::lag(ND))%>%plot_time_series(TIMESTAMP,diff_ND,.smooth=F)